Tensor network simulation of phase diagram of frustrated J1-J2 Heisenberg model 

on a checkerboard lattice 
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We use the recently developed tensor network algorithm based on infinite projected entangled 
pair states (iPEPS) to study the phase diagram of frustrated antiferromagnetic J1-J2 Heisenberg 
model on a checkerboard lattice. The simulation indicates a Neel ordered phase when J2 < O.88J1, a 
plaquette valence bond solid state when 0.88 < J2/J1 < 1.11, and a stripe phase when J2 > l.llJi, 
with two first-order transitions across the phase boundaries. The calculation shows the cross-dimer 
state proposed before is unlikely to be the ground state of the model, although such a state indeed 
arises as a metastable state in some parameter region. 
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Understanding frustrated quantum magnetic models 
is a long standing difficult problem in strongly correlated 
physics. Theoretical tools to study these systems are lim- 
ited. Exact diagonalization is limited by the small system 
size, and quantum Monte-Carlo simulation is hindered by 
the infamous sign problem. Among the frustrated mod- 
els, the antiferromagnetic J1-J2 Heisenberg model on a 
checkerboard lattice (or called the crossed chain model) 
is an important example that has raised a lot of interest 
[11-01 , due to its rich phase diagram and connection with 
real materials. This model is described by the Hamilto- 
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where Ji is the nearest-neighbor spin coupling rate on 
a square lattice and J2 is the next-nearest-neighbor spin 
coupling rate on a checkerboard pattern of plaquettes, 
as depicted in Fig. 1. A number of works provide com- 
plimentary studies in different parameter regions. The 
complete phase diagram for this model, however, still re- 
mains controversial. It is known that the system has 
a collinear long-range Neel order when J2 << Ji- At 
Ji « J2, the calculation based on the strong-coupling 
expansion predicts a plaquette valence bond solid as the 
ground state [H, 0-01 ■ In the region with J2 > Ji > 0, 
the phase is still under debate. Possibilities include the 
fourfold degenerate state with longrange spin order, sup- 
ported with semi-classical studies jGl] and large- expan- 
sion calculation Q, the sliding Luttinger liquid phase, 
supported with perturbative random phase calculation 
and exact diagonalization of a small system [ij , and the 
cross dimer state, supported with bosonization approach 
and two-step DMRG (density-matrix renormalization 
group) simulation [4]. 

Recently, tensor network algorithms emerge as a 
promising method to solve two-dimensional frustrated 
quantum systems [8-12[. There are different types of 



tensor network algorithms, but all the algorithms share 
the basic idea to describe the ground state of the model 
Hamiltonian as a tensor network state that respects the 
entanglement area law. The tensor network algorithms 
belong to the variational method and have no intrinsic 
sign problem for frustrated systems. The tensor net- 
work algorithms have been tested for a number of non- 
frustrated Hamiltonians, and the results agree pretty well 
with quantum Monte Carlo simulation [10|. Recently, 
the algorithms have also been applied to the frustrated 
Heisenberg model on a Kagome lattice ll| and the Ji- 
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J2-J3 model on a square lattice 

In this paper, we use a particular type of tensor net- 
work algorithm, the iPEPS (infinite PEPS) ji], to sim- 
ulate frustrated antiferromagnetic J1-J2 model on a a 
checkerboard lattice in the thermodynamic limit. We 
construct the complete phase diagram with the following 
findings: (1) the simulation shows two first-order phase 
transitions respectively at J2/J1 = 0.88 and J2/J1 = 
1.11, first from a Neel state to a plaquette valence bond 
solid, and then to a spin ordered stripe phase. (2) In the 
region with J2/J1 > 1.11 (except for the special point 
Ji = 0), our calculation supports the four-fold degener- 
ate states proposed in Ref. 6] as the ground state. In 
particular, the spin stripe phase seems to be the most 
stable one under perturbation. (3) Our simulation pro- 
vides strong evidence to show that the cross-dimer state 
is not the ground state of the system, although it indeed 
emerges as a metastable state in some parameter region 
(its energy is always significantly higher compared with 
the spin stripe phase). 

In the iPEPS algorithm, at every site, we represent the 
state as a five-index tensor with one physical index (with 
dimension two for a spin-half system) and four virtual 
indices (with internal dimension denoted by D) The 
wave function can be obtained by contracting all the vir- 
tual indices. To obtain the expectation value of a physical 
quantity, we need to first contract the physical index to 
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FIG. 1: Illustration of the Ji — J2 Heisenberg model on a 
checkerboard lattice. 

form a tensor network with internal dimension D^, and 
this tensor network is then contracted from the infinite 
boundary through multiplication with a matrix product 
state with internal dimension x- We tune the value of x 
until convergence is achieved in the calculated physical 
quantity. As a thumb of rule, typically at x ^ the 
relative error of energy induced by variation of x has been 
reduced to the order of 10^^, which indicates good con- 
vergence already. The dominant error of the calculation 
is from the small value of the internal dimension D. As 
the calculation time scales with D as D^^ (under choice 
of X ~ D'^) 0, the value of D in our simulation is limited 
to be about 5. We compare the energy as well as several 
other quantities (including the phase boundary specified 
below in Fig. 2) calculated with D — 4 and D — 5, and 
the difference is within a percent level. As an estimate, 
we expect that the relative error of our numerical sim- 
ulation, dominated by the limited value of D, is within 



or about a percent level for any short-range correlation 
function. 

The original iPEPS algorithm needs to assume transla- 
tional symmetry for calculation in the thermodynamical 
limit. The ground state of the Hamiltonian (1) can spon- 
taneously break the translational symmetry. To take into 
account the spontaneous symmetry breaking, we take a 
large unit cell and assume the translational symmetry 
only among different cells with no symmetry restriction 
for the tensors within the cell. In our simulation, the 
unit cell has 4x4 sites which is large enough to incor- 
porate the relevant ground states for this Hamiltonian 
that break the translational symmetry [l^. We apply 
imaginary time evolution to reach the ground state of 
the Hamiltonian. To avoid being stuck in a metastable 
state, we take a number of random initial states for the 
imaginary time evolution and pick up the ground state 
as the one which has the minimum energy over all the 
trials. 

In Fig. 3, we show the complete phase diagram for 
the Hamiltonian (1) from this calculation. To charac- 
terize the phase transition, we calculate derivative of the 

ground state energy = ^{{i j)) ('^' ' '^') ^^^^ respect 
to J2 (Ji is taken as the energy unit) and identify the 
singular point of this derivative as the phase transition 
point. To characterize properties of different phases, we 
calculate the spin order parameter ^5^^ for all sites i and 
the plaquette order parameter Qr.fl^A[14|. defined by 
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where a, /3, 7, and S denote four spins on a plaquette 
as shown in Fig. 1. Different phases are associated with 
different characteristic values of these parameters. For 
instance, a spin ordered state is characterized by a sig- 
nificant value of ^"^i^; in contrast, the plaquette valence 
bond solid state is characterized by a near-unity Qa^jS 
and a vanishing ( Si 



In the inserts of Fig. 2, we show the order parame- 
ters (^^i^ and QapyS as functions of J2/Ji- These order 
parameters change abruptly at the corresponding phase 
transition points, and the points of abrupt change are in 
agreement with the singularity points of The or- 

der parameters and the derivative of the ground state 
energy both have finite jumps at the phase transition 
points, which strongly indicates that we have two first- 



order transitions as we increase the ratio J2/ Ji, first from 
a Neel ordered state to a plaquette valence bond solid 
state at J2/J1 — 0.88, and then from the valence bond 
solid state to another spin-ordered phase (its nature will 
be discussed below) at J2/J1 = l-H. The possibility of 
two second order phase transitions with a coexistence re- 
gion of the spin and the valence bond solid orders in the 
intermediate region has been discussed in the literature 
[l5| . Within the resolution of our numerical simulation 
(the resolution is 0.01 for the ratio J2/J1 near the phase 
transition points), we do not find a coexistence region 
and the result supports a direct first-order transition. 

The nature of ground states in these three phases are 
further studied through calculation of the spin correlation 
function. In Fig. 3, We show the nearest-neighbor spin- 
spin correlation ( Si • Sj ) and orientation of local spins 
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FIG. 2: (Color online) The main plot shows dE/dJ2 as a 
function of J2/(Ji + J2)- Insets show the plaquette order 
(solid line with crosses) and the spin order (red dashed line 
with circles) as functions of J2/Jinear the transition points. 
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FIG. 3: (Color online) The upper left figure (a) shows orienta- 
tions of local spins (red arrows) on the checkerboard lattice at 
J2 = 0.5. The upper right figure (b) shows nearest neighbor 
spin-spin correlations < Si ■ Sj > along horizontal, vertical 
and diagonal bonds on a 16-site unit-cell at J2 — 1.0. The 
width and colors of the bonds are scaled such that the neg- 
ative correlation is represented by thicker bonds with darker 
color. The lower figures show both spin orientation and corre- 
lations < Si ■ Sj > in a Neel* state (c) and a stripe phase (d) 
at J2 = 2. Fig. (e) illustrates the dimer state, which appears 
as a metastable state at some parameters. 



with respect to the first spin on the up-left corner in each 
16-site unit-cell for three different J2 values. At J2 = 1, 
strong nearest-neighbor spin correlations (valence bonds) 
around the plaquettes breaking the lattice translational 
symmetry suggests a plaquette ordered state near this 
point, consistent with the finding in Fig. 2. At J2 — 0.5, 
the spin orientation indicates a conventional Neel ordered 
state. At J2 = 2, antiferromagnetic Neel order appears 
along the diagonal chains, but not in the horizontal or 
vertical axes. With imaginary time evolution starting 
from randomly chosen initial states, we actually get two 



different kinds of spin configurations shown in Fig. 3(c) 
and 3(d) for the final state. Their energies are almost 
degenerate within resolution of our numerical program. 
These spin configurations are in agreement with the four- 
fold degenerate states found in Ref. Q based on the 
large-iS* expansion (the other two degenerate states are 
obtained from Fig. 3(c) and 3(d) through a 90-degree 
rotation of the spin orientation). The spin configuration 
in Fig. 3(c) is called the Neel*-state in Ref. [3|, where 
the single-site spin '[~ or 4- in the conventional Neel state is 
replaced by the two-site unit tt or H. The configuration 
in Fig. 3(d) represents a spin stripe phase, where the spin 
orientations form a stripe along the horizontal or vertical 
direction, breaking the lattice rotational symmetry. The 
stripe phase is also predicted in the large- calculation 

a. 

To further clarify the phase at J2 = 2, we also show 
the spin-spin correlation (^Si ■ Sj^ in Fig. 3(c) and 3(d). 
Strong nearest-neighbor spin correlation appears along 
the diagonal chains. However, the distribution of these 
spin correlations does not break the symmetry of the lat- 
tice, so it is not a cross dimer or other valence bond solid 
state. The cross dimer state is predicted for this model in 
[1,13] for some region of J2I J\- In a cross dimer state, the 
nearest-neighbor spin correlations form the cross dimer 
pattern illustrated by Fig. 3(e). We indeed get this kind 
of cross dimer configuration from the imaginary time evo- 
lution starting from a pure cross dimer state for a certain 
region of J2/J1 as shown in Fig. 4. However, the ener- 
gies of the cross dimer states are strictly higher than the 
four-fold degenerate states shown in Fig. 3(c) and 3(d). 
So the cross dimer state is only a metastable phase in 
this region and does not give the real ground state. We 
compare the energy of our calculation with the energy 
of the DMRG calculation in Ref. [3], and our energy is 
significantly lower than the corresponding result in 
on the side J2 > Ji. For instance, at J2 = 2, our result 
shows a ground state energy oi E = —0.876 Ji for a stripe 
state, much lower than the energy of _B ~ — 0.75 Ji for a 
cross dimer state at the corresponding point in Ref. Q. 
Because of this large energy difference, it is unlikely that 
the cross dimer state emerges as the real ground state of 
the system. 

Some literature predicts a sliding Luttinger liquid state 
on the J2>Ji side of the antiferromagnetic checkerboard 
model [l|, 0] . We do not find evidence to support a tran- 
sition to a sliding Luttinger liquid state in our numerical 
simulation. Of course, due to the limitation of the inter- 
nal dimension D of the variational tensor network state 
in our calculation, it is possible that the sliding Luttinger 
liquid state is poorly approximated by the tensor network 
state with a small internal dimension and thus missed in 
our numerical simulation. We can not rule out this pos- 
sibility, however, we feel its chance is small due to the 
following test: we know at the limiting point Ji = 0, 
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FIG. 4: (Color online) Ground state energy calculated by 
iPEPS (blue squares) for different J-ij^J^ + J\)- The dashed 
line denotes energy of a pure dimer state, the green dia- 
monds denote energy of meta-stable dimer states calculated 
by iPEPS with imaginary time evolution from an initial pure 
dimer state. 
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FIG. 5: (Color online) (a) Long-range spin correlations along 
a diagonal chain obtained with D = 3 (crosses), D = 4 
(squares), and D = 5 (circles) at J2 = 2. (b) The same corre- 
lation for a Heisenberg chain calculated with D — 3 (circles), 
D = 5 (squares), D = 15 (diamonds), and D = 30 (triangles). 

the model reduces to decoupled Heisenberg chains whose 
ground state is described by a Luttinger liquid with al- 
gebraic decay of the spin correlation function. We use 
the same tensor network algorithm to calculate the long 
range spin correlation for the limiting case at Ji = 0. At 
this ID limiting point, we can have a much larger internal 
dimension D in numerical simulation, and in Fig. 5(b) 
we compare the result with D varying from 2 to 30. We 
see that the result at D = 5 has correctly demonstrated 
the algebraic decay of the spin correlation function and 
almost converged to the result a.t D — 30. So we do not 
necessarily need a large internal dimension for the tensor 
network algorithm to uncover the algebraic decay associ- 
ated with a spin liquid state. Keeping the same internal 
dimension at _D = 5, we turn on Ji (now a 2D model 
with J2 > Ji > 0), and find that long range spin corre- 
lation appears along the diagonal chains (see Fig. 5(a)), 
indicting that the spin order along this direction is very 
likely a real effect. 

Although our numerical program can not distinguish 
the energy of the four degenerate states shown in Fig. 
3(c,d) at the J2 > Ji side, it is very likely that the stripe 



phase will emerge as the real ground state in practice 
because of its robustness to perturbation in the Hamil- 
tonian. In real realization of the model Hamiltonian (1), 
the Ji coupling along the horizontal and the vertical di- 
rections might be slightly different, or apart from the 
J2 coupling on the checkerboard pattern, there might be 
small antiferromagnetic J2 coupling along the other pla- 
quettes. With any of these types of perturbation (which 
sound to be inevitable in practice), the energy of the 
Neel*-state will be lifted, and the stripe phase will emerge 
as the unique ground state of the system. 

In summary we have used the iPEPS method, a type 
of tensor networks algorithms, to calculate the ground 
states of the frustrated anti-ferromagnetic J1-J2 Heisen- 
berg model on a checkerboard lattice. We construct the 
complete phase diagrams, indicting two first order phase 
transitions, first from a Neel state to a plaquette valence 
bond solid and then to a spin stripe phase. The calcula- 
tion helps to clarify some of the previous debates on the 
phase diagram of this important model and provides a 
novel example for applications of the recently developed 
tensor network algorithms to frustrated systems. 
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